%cd '/Users/tessabold/Dropbox/Tessa and Tobi/Programmes'
clear all
close all
% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the standard limited commitment village
% economy for both general and partial equilibrium

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version June 2021
% =====================================================================

global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table AGG vs vss ScaleRov

filetosave='Vilcoal_replicate_June_2018';
write=1; % set this to 1 if you want to save tables, workspace and figures
%====================================================
% I. General equilibrium
% ====================================================
global P S gammagrid Y AGG
global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table AGG vs vss ScaleRov

nyvec=[3]; % number of indiv income values
rho=1.0;%linspace(,2,1); % laczo (2012, Wp) has 1 - log-preferences
 
  datapath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Data\Output';
  programpath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Programmes';

      relvarincscale=1;  % determines how to scale the difference in consumption growth variance
      analytic =1; % analytic VCV used - not needed
Pun_fac=0.0; % laczo (2012, Wp) has 0
datafile='armoments_fe0';
cd(datapath)
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);

cd(programpath)
COUNT=0;
village=2;

rho1=covyylag(village)/vary(village)

vareps=(vary(village)).*(1-rho1.^2);
if village==1
    N=34-1;                 % the number of total village members is N+1 under individual deviations
    beta=0.89;

elseif village==2;
    N=37-1;                 % the number of total village members is N+1 under individual deviations
      beta=0.90;
      N=180-1;
      betavec0=[0.9 0.94];
elseif village==3
    N=31-1;                 % the number of total village members is N+1 under individual deviations
      beta=0.88;
end
betavec=betavec0;
Nvill=1;
Nincproc=1;         varnu=0;                % number of income processes to look at
incproc=1;
T=1200; Tinit=201;                 % number of simulation periods in total, and number of initial periods to discard
maxcerr=0;                          % max cons meas error in log levels
N_cerr=2;                           % number of support points for cons meas error
maxincerr=2/3;                      % max inc meas error in multiples of inc variance (in levels)
TT=50000; % number of periods
vss=N;
nyagg=5; % number of agg inc values
% uninteresting options
ScaleRov=vss; % the scaling factor that transforms rov aggregate into individual income
Qrho=1; % number of rho parameters
Qbeta=1;% number of beta parameters
% set this to 1 if you want to calculate the moments for levels of
% consumption based on log consumption
momentclevinlog=1;


for ny=nyvec
COUNT=COUNT+1;
% 1. Income Process
% specify a markov process without persistence for individual income
% now we use Rouwenhorst
[incomevalsalt, Q1] = rouwenhorst(rho1,vareps^0.5,ny);
incomevals=reshape(exp(incomevalsalt'),length(incomevalsalt),1);

%  Simulate HH income for N HH
QQ1=Q1^1000;
SSdist=QQ1(1,:);
Qsum=cumsum(Q1);
for jj=1:N
    rng(jj,'twister');
    z(jj,1)=min(find(rand(1,1)<=cumsum(SSdist)));
    y_sim(jj,1)=incomevals(z(jj,1));
end
for n=1:N
    for t=2:TT
        rng(n*1000+t,'twister');
        z(jj,t)=min(find(rand(1,1)<=cumsum(Q1(z(jj,t-1),:))));
        y_sim(n,t)=incomevals(z(jj,t));
    end
end

% estimate aggregate income process
XX=[ones(1,TT-100)',log(sum(y_sim(:,100:end-1),1))'];
YY=log(sum(y_sim(:,101:end),1))';
beta_rov=inv(XX'*XX)*XX'*YY;
eps_Y=(YY'-beta_rov'*XX')';
Var_eps_Y=(eps_Y'*eps_Y)/TT;
[S_rov_temp, P_rov] =rouwenhorst(beta_rov(2),Var_eps_Y^0.5,nyagg);
S_rov=exp(beta_rov(1)/(1-beta_rov(2))+S_rov_temp)';
S_rov=S_rov/ScaleRov;

% ======== c. state of the economy: combinations of indiv and aggregate income ============
S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
P=kron(Q1,P_rov);
cumP=cumsum(P')';
cumQ=cumsum(Q1')';
state=length(S);
Y=S(:,1)+ScaleRov*S(:,2);			% possible aggregate income outcome
% ============ transition matrix for state of the economy ============
agent=2;
for Qbeta=1:length(betavec)
    beta=betavec(Qbeta);
clear caut eeewaut waut ewperf share normperf eewperf wsbnewsol wsb 
% Compute set of autarky values
for k=1:agent
    caut(:,k)=S(:,k)*(1-Pun_fac);
end
waut=zeros(state,agent);    

% autarky values
eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
% perfect insurance values
ewperf(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,Y./(vss+1));
share=linspace(0.3/(1+ScaleRov),2/(1+ScaleRov),100); % constant share of individual
for i=1:100
    ewperf(:,1,i)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,Y*share(i));
    ewperf(:,2,i)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,Y*(1-share(i))/ScaleRov);
    normperf(i)=max(max(abs((eeewaut(:,:)-ewperf(:,:,i)))));
end
% check if perfect insurance is feasible
[mintemp,minind]=min(normperf);
eewperf=ewperf(:,:,minind);
clear ewperf mintemp minind
if eewperf-eeewaut>0
    disp('perfect insurance')
  %  pause
end

for j=1:size(S,1)
    waut(j,1)=utility(rho,caut(j,1))+beta*P(j,:)*eeewaut(:,1);
    % autarky values for the rest of village are always the same - perfect
    % risk-sharing among the rest of village
    waut(j,2)=utility(rho,caut(j,2))+beta*P(j,:)*eeewaut(:,2);
    
end

% ====================================================
% II. Compute constrained insurance
% ====================================================

% grid of time-varying planner weights
lambda=[0.005,0.1,0.2,0.25,0.3,0.35:0.05:0.65,0.7,0.75,0.8,0.9,0.995]';
m=length(lambda);
gammagrid=[lambda vss/ScaleRov*(1-lambda)];

gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];

% ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
n=length(gammagrid);
w=zeros(n,state,agent);
gammar=zeros(n,state,agent);
c=zeros(n,state,agent);
phi=zeros(n,state,agent);

% consumption given aggregate resources and planner weights
for j=1:state
    % for k=1:agent
    c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
    c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
    % end
end
cfirstbest=c;
%  relative planner weights
for j=1:state
    for k=1:agent
        gammar(:,j,k)=gammagrid(:,agent+k);
    end
end

% Compute the first best value functions without binding
% constraints
w=zeros(n,state,agent);
%T first best expected utility
for j=1:state
    w1(:,j,1)=utility(rho,cfirstbest(:,j,1))+beta/(1-beta)*[utility(rho,cfirstbest(:,:,1))]*P(j,:)';
    w1(:,j,2)=utility(rho,cfirstbest(:,j,2))+beta/(1-beta)*[utility(rho,cfirstbest(:,:,2))]*P(j,:)';
end
dww=1;
wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
while dww>0.00001
    w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
    w(:,:,2)=utility(rho,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
    
    dww=max(max(max(abs(w-wold))));
    wold=w;
    
end
wfirstbest=w;

% continuation utility as function of planner weight and state of the
% economy
for k=1:m
    for j=1:state
        ewfirstbest(k,j,1)=P(j,:)*wfirstbest(k,:,1)';
        ewfirstbest(k,j,2)=P(j,:)*wfirstbest(k,:,2)';
    end
end
for k=1:m
    for j=1:state
        ewautaut(k,j,1)=P(j,:)*waut(:,1);
        ewautaut(k,j,2)=P(j,:)*waut(:,2);
        wautaut(k,j,1)=waut(j,1);
        wautaut(k,j,2)=waut(j,2);
    end
end

% ====== a. Start the loop with the enforcement constraints ======
%%%this will be the parameterized expectation
ewsb=ewfirstbest;
wsb=wfirstbest;
gap=1;
options=optimset('Display','off');
count=0;
eps=0.0000001;
t=0;
ggammagrid=fix(gammagrid*100);
itgap=0;
maxitgap=200;
while gap >=0.01 && itgap<=maxitgap
    clear rr
    itgap=itgap+1;
    coalnonsust=zeros(n,state);
    count=0;
    % Begin proper loop
    
    count=0;
    countgrid=[ones(n,state)];
    for j=1:state
        %%no permutations here
        disp(j)
        for i=1:n
            count=count+1;
            index=[i,j];
            
            % Now start updating all the functions c, z, w
            
            % 1. Constraint binding for Person 1 only
            if utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)<waut(j,1) && utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)>=waut(j,2) ;
                bind=1; %%i.e. agent 1 has the binding constraint
                rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1));
                if isempty(rr)==1
                    coalnonsust(i,j)=1;
                     nonconv(i,j)=nan;
                else
                    %T x0 is the guess for consumption of constrained agent
                    %1, corresponding to planner weight rr(1)
                    x0(1) = cfirstbest(rr(1),j,1);
                    %T vectors of consumption and utiltiy that meet
                    %participation constraints
                    
                    [x, output]=agent1maxrho_rov_15_Mar(x0,index,ewsb,waut,bind);
                    
                    %T 29_Mar This keeps track of nonconvergence, but
                    %allows the programme to continue - sometimes it
                    %converges in a subsequent iteration
                    if output(3) ~=1
                        nonconv(i,j)=output(3);
                        wsbnewsol(i,j,1:agent)=wsb(i,j,1:agent);
                    else
                        nonconv(i,j)=nan;
                        csol(i,j,1)=x(1);
                        gammar(i,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                        csol(i,j,2)=(Y(j)-x(1))/ScaleRov;
                        wsbnewsol(i,j,1:agent)=output(1:agent);
                        phisol(i,j,1)=gammagrid(i,4)^rho/gammar(i,j,2)^rho-1;
                        phisol(i,j,2:agent)=0	;
                    end
                    
                    
                     
                end
                
                % 2. Constraint binding for person 2 only
            elseif utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)>=waut(j,1) && utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)<waut(j,2);
                bind=2; %%i.e. agent 2 has the binding constraint
                rr=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                if isempty(rr)==1
                    coalnonsust(i,j)=1;
                     nonconv(i,j)=nan;
                else
                    x0(1) = cfirstbest(rr(end),j,1);
                    [x, output]=agent1maxrho_rov_15_Mar(x0,index,ewsb,waut,bind);
                    
                    
                    if output(3)~=1
                        nonconv(i,j)=output(3);
                        wsbnewsol(i,j,1:agent)=wsb(i,j,1:agent);
                        
                    else
                        nonconv(i,j)=nan;
                        csol(i,j,1)=x(1);
                        gammar(i,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                        csol(i,j,2)=(Y(j)-x(1))/ScaleRov;
                        wsbnewsol(i,j,1:agent)=output(1:agent);
                        
                        phisol(i,j,2)=gammar(i,j,2)^rho/gammagrid(i,4)^rho-1;
                        phisol(i,j,1)=0;
                    end
                    
                    
                end
                
                % 3. Constraint binding for none
            elseif utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)>=waut(j,1) && utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)>=waut(j,2) ;
                csol(i,j,1:agent)=cfirstbest(i,j,1:agent);
                gammar(i,j,1:agent)=gammagrid(i,agent+1:2*agent);
                phisol(i,j,1:agent)=0;
                
                for g=1:agent
                    wint(g)=interp1(gammagrid(:,4),ewsb(:,j,g),gammagrid(i,4));
                end;
                wsbnewsol(i,j,1)=utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1);
                wsbnewsol(i,j,2)=utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2);
                nonconv(i,j)=nan;
            end % of if loop
        end  % end of i loop
        
    end  % end of j loop
    
    if sum(sum(coalnonsust))>0 || sum(sum(isnan(nonconv)))<state*n  && itgap>=3
        % break
        % disp('break')
        coalnonsustind(Qrho,Qbeta,vss)=mean(mean(coalnonsust(find(coalnonsust>0))));
        nonconvind1(Qrho,Qbeta,vss)=sum(sum(isnan(nonconv)));
        gap=100;
    else
        %T 29_Mar is the norm taken over agent 1 only?
        gap=norm(wsbnewsol(:,:,1)-wsb(:,:,1))
    end
    nonconvind=mean(nonconv(find(~isnan(nonconv))));
    if min(min(nonconv))<0
        disp('Problem with convergence in this iteration')
        %pause
        %Nonconvind(vss,Qrho,Qbeta)=mean(nonconv(find(nonconv(:,:)>0,vss)))
    end
    
    wsb=wsbnewsol;
    
    for k=1:m
        for j=1:state
            ewsb(k,j,1)=P(j,:)*wsb(k,:,1)';
            ewsb(k,j,2)=P(j,:)*wsb(k,:,2)';
        end
    end
    
end % of iteration loop

% =========================================================
% Simulate the economy
% =========================================================
clear csim Yagg income indSjj incomeind gammarfun phinew
                csim=nan(vss,T);
                Yagg =nan(1,T);
                income =nan(vss+1,T);
                t=1;
                %initial state
                rng(3,'twister');
                indSjj(:,1)=randi(ny,N+1,1);
                   
                   
                    for k=1:vss+1
                        rng(k,'twister');
                            incomeind(k,t)=randi(ny);
                            income(k,t)=incomevals(incomeind(k,t));  
                            csim(k,t)=income(k,t);
                    end
                        Yagg(t)=sum(income(:,t));
                        [Smintemp(t),Sind(t)]=min(abs(Yagg(t)/(vss+1)-S(1:nyagg,2)));
                        for k=1:vss+1
                           indSjj(k,t)= (incomeind(k,t)-1)*nyagg+Sind(t);
                        end
                        
                while t<T
                     t=t+1
                      for k=1:vss+1
                          rng(t*100+k,'twister');
                            incomeind(k,t)=min(find(cumQ(incomeind(k,t-1),:)>=rand(1,1)));
                            income(k,t)=incomevals(incomeind(k,t));  
                    end
                        Yagg(t)=sum(income(:,t));
                        [Smintemp(t),Sind(t)]=min(abs(Yagg(t)/(vss+1)-S(1:nyagg,2)));
                        for k=1:vss+1
                           indSjj(k,t)= (incomeind(k,t)-1)*nyagg+Sind(t);
                        end
                   

                    for k=1:vss+1
                        
                        
                        % check that aggregate and individual incomes are
                        % consistent with the states caculated above
                      
                       % calculate weights at beginning of period
                        gammarfun(k,t)=(Yagg(1,t-1)-csim(k,t-1))/csim(k,t-1)/ScaleRov;
                        gammarfun(k,t)=min(gammarfun(k,t),max(gammagrid(:,4)));
                        % calculate the new planner weight for all individuals
                        phinew(k,t)=interp1(gammagrid(:,4),phisol(:,indSjj(k,t),1),gammarfun(k,t));
                    end

                    phinew(:,t)';
                    x0=csim(:,t-1)*Yagg(1,t)/Yagg(1,t-1);
                    vs=vss+1;
                    % solve for consumption
                    [x,output]=conssolve(x0,phinew(:,t),csim(:,t-1),Yagg(1,t));
                    csim(:,t)=x';
                    if isnan(csim(:,t))
                        stop
                    end
                    
             
end


% =========================================================
% calculate moments of consumption and income growth
% =========================================================
clear incgrowth cgrowth cshare_raw dcshare_raw 
csim0=csim;income0=income;
csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
            rng(5,'twister'); cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
            rng(6,'twister'); incerr=randn(size(csim0));   Momentsest1=nan(30,1); ind=[];
            % loop over size of cons measurement error
           for kkk=1:size(csim0,1)
                tempvar(kkk)=var(log(csim0(kkk,:)));
            end
            Tempvar=mean(tempvar);
          ic=1;
                csim=exp(log(csim0));
                logcsim=log(csim);
                logcsim0=log(csim0);
                logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
                logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
                % add inc meas error back in
                income=exp(log(income0)+varnu(incproc)^0.5*incerr);
                logincome=log(income);
                logincome0=log(income0);
                logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
                income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
                
                
                for ivill=1:Nvill
                    aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    logaggvillincome=log(aggvillincome);
                    
                    logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    
                    
                end
                aggincome=sum(income,1);
                
                
                
                % === discard first Tinit-1 periods
                % this chooses log or level consumption for
                % moments
                if momentclevinlog==1
                    ccc=logcsim(:,Tinit:T);
                    ccc_cond=logcsim_meandev(:,Tinit:T);
                    ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
                else
                    ccc=(csim(:,Tinit:T));
                    ccc_cond=(csim_meandev(:,Tinit:T)); 
                end
                yyy=logincome(:,Tinit:T);
                yyy_cond=logincome_meandev(:,Tinit:T);
                yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
                
                cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
                cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
                cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
                
                yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
                yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
                aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
                
                
                vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
                
                % calculate moments for each individiual
                for ind=1:size(csim,1)
                    
%                     dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
%                     dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
%                     dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
%                     dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
%                     
%                     % conditioning on income increases, and income
%                     % decreases
%                     % NB: we allocate no change to decrease!!
%                     rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
%                     rrhc=find(yy(ind,:)>0);
%                     rrlc=find(yy(ind,:)<=0);
%                     rrhc_cond=find(yy_cond(ind,:)>0);
%                     rrlc_cond=find(yy_cond(ind,:)<=0);
%                     rrhc_groupcond=find(yy_groupcond(ind,:)>0);
%                     rrlc_groupcond=find(yy_groupcond(ind,:)<=0);
%                     
%                     
%                     %%%THE VARIANCES
%                     
%                     vardclog(ind,vss)=var(cc(ind,:));
%                     vardylog(ind,vss)=var(yy(ind,:));
%                     vardc_dypos(ind,vss)=var(cc(ind,rrhc));
%                     vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
%                     
%                     %%%conditional on village income!
%                     vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
%                     vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
%                     vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
%                     vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));
%                     
%                     %%%conditional on group income!
%                     %Step 1: create your residuals
%                     vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
%                     vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
%                     vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
%                     vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
%                     vardc_dypos(ind,vss)=var(cc(ind,rrhc));
%                     vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
%                     
%                     %%%and the relative variances
%                     
%                     reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));
%                     
%                     reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
%                     reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);
%                     
%                     
%                     %%%THE BETA COEFFICIENTS
%                     
%                     %1) Unconditional
%                     %risk-sharing whole sample
%                     covtemp=cov(cc(ind,:),aggyy(1,:));
%                     betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
%                     covtemp=cov(cc(ind,:),yy(ind,:));
%                     betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
%                     %risk-sharing for those with income gains versus
%                     %those with income losses
%                     YY=cc(ind,:)';
%                     yy_hc(ind,:)=zeros(1,size(cc,2));
%                     yy_hc(ind,rrhc)=yy(ind,rrhc);
%                     oneone=ones(1,size(cc,2));
%                     ones_hc(ind,:)=zeros(1,size(cc,2));
%                     ones_hc(ind,rrhc)=oneone(rrhc);
%                     XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
%                     beta_yhigh(:,ind,vss)=regress(YY,XX);
%                     betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
%                     betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
%                     
%                     %2) conditional on villag income
%                     covtemp=[];
%                     X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
%                     YY=cc_cond(ind,:)';
%                     beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);
%                     
%                     X=[ones(size(aggyy(1,:),2),1),aggyy'];%
%                     YY=cc(ind,:)';
%                     beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);
%                     
%                     %risk-sharing for those with income gains versus
%                     %those with income losses
%                     YY_cond=cc_cond(ind,:)';
%                     yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
%                     yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
%                     oneone=ones(1,size(cc_cond,2));
%                     ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
%                     ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
%                     XX_cond=[ones(1,size(cc_cond,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
%                     beta_yhigh_cond(:,ind,vss)=regress(YY_cond,XX_cond);
%                     
%                     %%
%                     
%                     %2) conditional on group income
%                     covtemp=[];
%                     X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
%                     YY=cc_groupcond(ind,:)';
%                     beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);
%                     
%                     X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
%                     YY=cc(ind,:)';
%                     beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);
%                     
%                     %risk-sharing for those with income gains versus
%                     %those with income losses
%                     YY_groupcond=cc_groupcond(ind,:)';
%                     yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
%                     yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
%                     oneone=ones(1,size(cc_groupcond,2));
%                     ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
%                     ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
%                     XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
%                     beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);
                end 
                    XXreg=reshape(aggvillyy,size(aggvillyy,1)*size(aggvillyy,2),1);
                     XXreg=[ones(length(XXreg),1) XXreg];
                     YYreg=reshape(cc,size(cc,1)*size(cc,2),1);
                     bb=regress(YYreg,XXreg);
                     res_cc_cond=YYreg-XXreg*bb;
                     res_cc_cond=reshape(res_cc_cond,size(cc,1),size(cc,2));
                    
                    
                    
                    
                    
                    
                    

            RES_COND(1:size(res_cc_cond,1),:,ic,1,Qbeta)=res_cc_cond;
            DC_LOG_COND(1:size(cc_cond,1),1:size(cc_cond,2),ic,1,Qbeta)=cc_cond;
            DC_LOG(1:size(cc,1),1:size(cc,2),ic,1,Qbeta)=cc;
         


         
end
DC_LOG_COND2=DC_LOG_COND(:,:,1,:,:);
DC_LOG_COND2=squeeze(DC_LOG_COND2);
RES_COND2=RES_COND(:,:,1,:,:);
RES_COND2=squeeze(RES_COND2);
DC_LOG_2=DC_LOG(:,:,1,:,:);
DC_LOG_2=squeeze(DC_LOG_2);
save('permmodel2_village_RS.mat','DC_LOG_COND2','RES_COND2','DC_LOG_2')
end

 